Demonstrating how to run a multi-species Bayesian hierarchical model of bat fatalities at wind facilities across the United States

USGS Disclaimer

This software has been approved for release by the U.S. Geological Survey (USGS). Although the software has been subjected to rigorous review, the USGS reserves the right to update the software as needed pursuant to further analysis and review. No warranty, expressed or implied, is made by the USGS or the U.S. Government as to the functionality of the software and related material nor shall the fact of release constitute any such warranty. Furthermore, the software is released on condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from its authorized or unauthorized use.

Purpose

This vignette is intended to demonstrate how to run a multi-species Bayesian hierarchical model, developed at the U.S. Geological Survey, which predicts bat fatalities at U.S. wind facilities. We show how to load included data objects, run the model, and generate bat fatality predictions.

Introduction to the NationalBatsAndWind R Package

This repository is installed as an R package and has been primarily developed to support efficient data access and scientific reproducibility.

Data frames

Several prepared data frames are included in the package in order to easily run this analysis.

The first data frame describes turbines recorded in the U.S. Wind Turbine Database (USGS), including both those currently in operation and those that have been decommissioned. Load the data frame and its documentation using the following commands. The code to prepare these data can be accessed at: data-raw/wind/turbines.R file.

#Load the turbines data frame
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::turbines

#Run the following command to view the metadata for each column in `turbines`
#This will appear within the "Help" tab to the right in RStudio
?NationalBatsAndWind::turbines

The second data frame aggregates turbines by wind energy facilities. The data frame includes descriptive statistics for turbine specifications and other covariates associated with each facility. Facilities are grouped by development phase and decommission status. For example, if a facility has multiple phases, the facility may appear in multiple rows, but the start_year and end_year values for these phases will vary. This occurs because a single facility may have constructed and/or decommissioned wind turbines across several years. Load this data frame and its documentation using the following commands. The code to prepare these data are in the data-raw/wind/facilities.R file.

#Load the facilities data frame
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::facilities

#Run the following command to view the metadata for each column in `facilities`
?NationalBatsAndWind::facilities

Spatial grid

We processed spatial covariates on a 11.5 x 11.5 km grid across the contiguous U.S. The original grid was created by the National Renewable Energy Laboratory (NREL) to estimate and display wind energy development potential under a variety of conservation, social, and geographic restrictions. Processed covariate data across the grid is stored in data/grids/nrel/final_nrel_grid_with_covariates.shp. Use the following commands to load the processed grid and its documentation.

#Load the nrel_grid data frame
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::nrel_grid

#Run the following command to view the metadata for each column in `nrel_grid`
?NationalBatsAndWind::nrel_grid

To reduce the package size, the raw data layers to process these covariates are not directly included in this repository. The code to reproduce this spatial grid can be found at data-raw/grid/add_covariates_to_nrel_grid.R. To rerun this script fully, please reach out to the point of contact for this repository to request copies of the raw data layers.

As an example of the data stored in the nrel_grid, we will plot average windspeed at 10 meters above ground across the continental U.S below.

#Load state data for outlines in map
state_lookup <- data.frame(state = state.abb, state_name = state.name)
states <- map_data("state")
states <- st_as_sf(states, coords = c('long', 'lat'), crs = 4326)
states <- st_transform(states, crs = st_crs(NationalBatsAndWind::nrel_grid))
states$x <- st_coordinates(states)[,1]
states$y <- st_coordinates(states)[,2]
states$state <- states$region

#Plot a continuous covariate from the NREL grid data
covariate_plot <- ggplot() +
    #Select windspeed at 10 m above ground
    geom_sf(
      data = NationalBatsAndWind::nrel_grid,
      aes(fill = windspeed_10_m_av),
      color = "transparent")  +
    #Add state outlines
    geom_polygon(
      data = states,
      aes(x = x, y = y, group = group),
      color = "white", fill = "transparent", linewidth = 0.25) +
    #Set graphical parameters for a pretty plot
    scale_fill_gradientn(
      colors = gist_ncar(100),
      labels = scales::unit_format(unit = "m/s", big.mark = ","),
      na.value="gray50") +
    ggthemes::theme_map() +
    labs(fill = "Average windspeed") +
    theme_minimal() +
    theme(plot.background = element_rect(fill = "white", color = "white"),
          legend.position = "right",
          text = element_text(size=16, colour="black"),
          axis.text=element_text(size=12, colour="black"),
          axis.title=element_text(size=12,face="bold")) +
    ylab("") +
    xlab("") +
    coord_sf(crs = crs(NationalBatsAndWind::nrel_grid))

#Save plot to folder
ggsave(
  "plots/nrel_grid_covariates_example.bmp",
  plot = covariate_plot,
  width = 9,
  height = 4,
  units = c("in"),
  dpi = 300,
  create.dir = T
)

#Print plot to vignette
covariate_plot

Fatality data

The bat fatality data in this software release originated from post-construction monitoring of bat fatalities at wind facilities across the United States. Because these data have been collected and delivered to the U.S. Fish & Wildlife Service (USFWS) related to the issuance of incidental take permits, there is often a requirement the data and results are available to the public. Each of the reports included in this release has previously been made publicly available, and we provide the public hyperlink for each report. The data in this release are formatted to support a reproducible analysis. Many other post-construction monitoring reports exist but are not included in this release due to business confidential rights/restrictions which may make them exempt from disclosure under Freedom of Information Act (FOIA) Exemption 4.

In our full analysis, which is not included here, we are able to run our model across all reports shared with the USFWS and USGS, but in this repository, we only include publicly available data. We plan to publish a manuscript on an analysis of our full data set soon. For now, we are releasing the software and a demo analysis on the subset of public data. These public data are stored in the public_pcm data frame and can be used to run a demo version of our model below. By showing how these data are formatted, we hope others with additional PCM data that we may not have access to will be able to prepare their PCM data in the format shown in public_pcm and run our model on their own data to independently verify our results. The scripts in data-raw can help prepare additional PCM data in the public_pcm format and here we relate covariate values from the spatial grid to the PCM reports.

#Load mortality data
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::public_pcm

#Check out the documentation for `public_pcm`
?NationalBatsAndWind::public_pcm

Model details

The following describes the model implemented in inst/jags/multispecies_nmixture_without_ecoregion_poisson.txt. The second version (inst/jags/multispecies_nmixture_with_ecoregion_poisson.txt) includes a random effect for ecoregion in the model. The Poisson-Binomial mixture model for the count of carcasses \(Y_{fs}\) observed at facility \(f\) for species \(s\), out of \(N_{fs}\) total estimated fatalities, is given by:

\[ \begin{aligned} Y_{fs} \ | \ N_{fs}, p_{f} &\overset{indep.}{\sim} \text{Binomial}\left(N_{fs}, p_{f}\right), \\ N_{fs} \ | \ \lambda_{fs} &\overset{indep.}{\sim} \text{Poisson}\left(\lambda_{fs}\right). \end{aligned} \]

A prior for the detection probability \(p_{f}\) is informed by the number of fatalities observed over the total bat fatalities originally estimated for facility \(f\) in its associated post-construction monitoring (PCM) report. Most of these PCM reports estimated total bats killed per megawatt, so we adjusted the estimates by the total rated capacity of the facility when applicable. These reports also used a variety of fatality rate estimators, some of which are known to be biased. Instead of using these potentially biased estimates as a fully trusted value, we borrowed the information to weakly inform the prior of the detection probability \(p_f\) at each facility. Parameters \(a_{f}\) and \(b_{f}\) were calculated to center the empirical detection probability of number of observed fatalities out of the total expected fatalities between the upper and lower uncertainty bounds calculated by the original estimator such that:

\[ p_{f} \sim \text{Beta}(a_{f}, b_{f}) \]

Then, based on the number of fatalities observed of each species at each facility and the weakly informed prior on the detection probability, the rate of fatalities \(\lambda_{fs}\) is modeled by a Bayesian GLMM with an offset based on the number of turbines at the facility \(n_{turbines}\):

\[ \begin{aligned} \log(\lambda_{fs}) = \text{log}(n_{turbines}) + \beta_s +& \\ & X_{fs,\text{seasons}}\beta_{\text{seasons}} + \\ & (1-I(\text{yes (1)/no (0) estimator adjusts for distance})) \times \frac{1}{\text{distance}_{f}} \beta_{\text{distance}} + \\ & X_{f, \text{covariates}}\beta_{\text{covariates,s}} + \\ & \epsilon_{\text{ecoregion}} + \epsilon_{\text{obs}} \end{aligned} \]

Each species has a species-specific fixed-effect intercept. Recall that JAGS parameterizes Normal distributions using mean and precision (a precision of 0.01 results in a very flat Normal distribution).

\[ \begin{aligned} \beta_{s} \sim \text{Normal}(0, 0.01) \end{aligned} \]

Each report has an adjustment for how much of the spring, summer, and fall seasons were searched. We assumed that if more of each season was searched, then more bat fatalities would be observed, so each \(\beta_{\text{seasons}}\) coefficient has a half-Normal prior with shared hyper-means and precision across species.

\[ \begin{aligned} \beta_{\text{season}_i, s} \sim \text{Normal}(\mu_{\text{season}_i}, \tau_{\text{season}}) T(0,) \end{aligned} \]

Some reports used an estimator that did not adjust for the number of bat carcasses that may have fallen outside of the search radius. For those reports, an adjustment is included in the model to account for the density of carcasses that may have fallen outside the searched area.

\[ \begin{aligned} (1-I(\text{yes (1)/no (0) estimator adjusts for distance})) \times \frac{1}{\text{distance}_{f}} \beta_{\text{distance}} \\ \end{aligned} \]

Species-specific coefficients for the covariate relationships \(X_{f, \text{covariates}}\beta_{\text{covariates,s}}\) were assumed to have covariate-specific hyper-means and precision.

\[ \begin{aligned} \beta_{\text{covariate}_j, s} \sim \text{Normal}(\mu_{\text{covariate}_j}, \tau_{\text{covariate}_j}) \end{aligned} \]

Additionally, there are optional random effects for ecoregion \(\epsilon_{\text{ecoregion}}\) and observation-level random effects \(\epsilon_{\text{obs}}\) to help account for overdispersion.

Other parameters in the model were assumed to have the following priors:

\[ \begin{aligned} \beta_{\text{distance}} &\sim \text{Normal}(0, 0.01) \\ \mu_{\text{season}_i} &\sim \text{Normal}(0, 0.01) T(0,) \\ \tau_{\text{season}} &\sim 1 / \sqrt{\text{Uniform}(0, 10)} \\ \mu_{\text{covariate}_j} &\sim \text{Normal}(0, 0.01) \\ \tau_{\text{covariate}_j} &\sim 1 / \sqrt{\text{Uniform}(0, 10)} \end{aligned} \]

Process input data

There are a few additional steps for processing the data that are required before running the model. The steps that we will process in the next steps include:

1. Determine whether facilities in the PCM data fall within each species range
2. Create a matrix and vector of relevant covariates
3. Center and scale the covariate matrix
4. Prepare the model matrix and other model inputs

Step 1: Species ranges

In this step, we join the species ranges with each PCM study to determine whether the facility falls within each species range. This step will help the model estimate species-specific covariate effects only when a facility falls within the expected range of the species or fatalities of the species were observed even though the facility is out of the known range.

#Load shapefile of species ranges
spc_ranges <- read_sf("../data-raw/bats/ranges/USGS_ScienceBase_coded.shp") %>% 
  st_transform(crs = st_crs(NationalBatsAndWind::nrel_grid))

#Determine species codes that exist in the public_pcm data
spc_start <- which(colnames(public_pcm) == "EPFU (Big brown)")
spc_end <- which(colnames(public_pcm) == "NYFE (pocketed fre-tailed)")
spcs <- colnames(public_pcm)[spc_start:spc_end]
spcs_code <- sapply(spcs, substr, 1, 4)

#Select the species ranges that exist within the public_pcm data
spc_ranges <- spc_ranges %>%
  filter(Spp_code %in% spcs_code) %>%
  st_transform(crs = st_crs(NationalBatsAndWind::nrel_grid))

#Add species ranges to the mortality data (1 if facility in range, 0 otherwise)
public_pcm_with_spc <- NationalBatsAndWind::public_pcm |>
  sf::st_as_sf(coords = c("x", "y"), crs = st_crs(NationalBatsAndWind::nrel_grid)) |>
  sf::st_join(spc_ranges) |> 
  dplyr::select(-SCI_NAME, -COMMON_NAM) |>
  tidyr::pivot_wider(
    names_from = Spp_code,
    values_from = Spp_code,
    values_fn = ~1,
    values_fill = 0) |>
  st_drop_geometry() |>
  #Join nrel_grid with the facility data
  dplyr::left_join(NationalBatsAndWind::nrel_grid, by = "fid")

Step 2: Covariates

Now we prepare the covariate data matrix for the facilities with PCM reports. In summary, a few of the covariates in the public_pcm are slightly altered or renamed to be consistent with the prediction steps later on. Additionally, several other covariate vectors that help with the bias adjustments for search strategies are created for the model.

#Create the covariate matrix
Xdf <- public_pcm_with_spc |>
  st_drop_geometry() |>
  #Edit values of a few covariates to account for zeroes and missing values
  dplyr::mutate(
    years_since_wns = uswtdb_min_p_year - wns_yoa_av,
    years_since_wns_nonneg = if_else(years_since_wns < 0, 0, years_since_wns),
    max_cutin_speed = suppressWarnings(as.numeric(max_cutin_speed)),
    adj_distance = grepl("no", estimator_accounts_for_max_dist) * 1,
    max_cutin_speed = if_else(is.na(max_cutin_speed), 3, max_cutin_speed),
    dist_searched_transform = 1/max_distance_searched) |>
  #Select covariates to include in model matrix
  dplyr::select(
    #Select covariates that are not scaled first
    fid, ecoregion2_mode, used_cutin,
    adj_distance, dist_searched_transform,
    prop_spring, prop_summer, prop_fall,
    n_turbine = estimate_relevant_tnum,
    #Rename and select these covariates
    total_rsa = estimate_relevant_rsa,
    total_mw = estimate_relevant_cap,
    #Also select the following
    starts_with("nlcd_ilr"),
    mean_rsa, mean_lrs, age,
    lrs_mean_windspeed_av, physdiv_av,
    vapor_pr_av
  )

#Some estimators account for fatalities that may fall outside of the search area
#This is a boolean vector that indicates which reports still need adjusting
adj_distance <- Xdf$adj_distance

Step 3: Center and scale covariate matrix

Here, we center and scale most of the covariates and save the parameter values for centering and scaling to use again during prediction.

#Center and scale the covariate matrix
Xdfcs <- Xdf |>
  #Center and scale all variables except those that are factors or indicators
  mutate(across(-c(1:9), \(x) scale(x, center = TRUE, scale = sd(x, na.rm = TRUE)))) |>
  #The season proportions don't handle scaling well. So center to 1 so that
  #a full season search drops out of the model (this is what we want the
  #model to estimate by default) and "scale" by 1 so no-scaling occurs
  mutate(across(starts_with("prop"), \(x) scale(x, center = 1, scale = 1))) |>
  #The distance_searched_transform has a half-normal prior, so let's set
  #the center to 0 (since estimates that don't need to be distance adjusted
  #will have a value of 0 and drop out of the model anyway) and scale
  mutate(across(dist_searched_transform, \(x) scale(x, center = 0, scale = sd(x, na.rm = TRUE)))) |>
  #Calculate quadratic effects for some variables after center/scale operation
  mutate(vapor_pr_av_sq = vapor_pr_av^2)

#Save the center and scale information for the covariate matrix to use later
cs_save <- data.frame(rbind(
  c('used_cutin', 0, 1),
  unlist(c('dist_searched_transform', attributes(Xdfcs$dist_searched_transform)[2:3])),
  unlist(c('prop_spring', attributes(Xdfcs$prop_spring)[2:3])),
  unlist(c('prop_summer', attributes(Xdfcs$prop_summer)[2:3])),
  unlist(c('prop_fall', attributes(Xdfcs$prop_fall)[2:3])),
  unlist(c('total_rsa',attributes(Xdfcs$total_rsa)[2:3])),
  unlist(c('mean_rsa',attributes(Xdfcs$mean_rsa)[2:3])),
  unlist(c('mean_lrs',attributes(Xdfcs$mean_lrs)[2:3])),
  unlist(c('age',attributes(Xdfcs$age)[2:3])),
  unlist(c('lrs_mean_windspeed_av', attributes(Xdfcs$lrs_mean_windspeed_av)[2:3])),
  unlist(c('physdiv_av',attributes(Xdfcs$physdiv_av)[2:3])),
  unlist(c('vapor_pr_av',attributes(Xdfcs$vapor_pr_av)[2:3])),
  unlist(c('vapor_pr_av_sq',attributes(Xdfcs$vapor_pr_av)[2:3])),
  unlist(c('nlcd_ilr_1',attributes(Xdfcs$nlcd_ilr_1)[2:3])),
  unlist(c('nlcd_ilr_2',attributes(Xdfcs$nlcd_ilr_2)[2:3])),
  unlist(c('nlcd_ilr_3',attributes(Xdfcs$nlcd_ilr_3)[2:3])),
  unlist(c('nlcd_ilr_4',attributes(Xdfcs$nlcd_ilr_4)[2:3])),
  unlist(c('nlcd_ilr_5',attributes(Xdfcs$nlcd_ilr_5)[2:3])),
  unlist(c('nlcd_ilr_6',attributes(Xdfcs$nlcd_ilr_6)[2:3])),
  unlist(c('nlcd_ilr_7',attributes(Xdfcs$nlcd_ilr_7)[2:3])),
  unlist(c('nlcd_ilr_8',attributes(Xdfcs$nlcd_ilr_8)[2:3])),
  unlist(c('nlcd_ilr_9',attributes(Xdfcs$nlcd_ilr_9)[2:3])))) |>
  mutate(
    scaled.center = as.numeric(scaled.center),
    scaled.scale = as.numeric(scaled.scale)) |>
  rename(colname = V1)

#Update rownames of the saved data and remove the names in the first column
rownames(cs_save) <- cs_save$colname
cs_save <- cs_save[,2:3]

#Some reports require adjustments for total amount of each season searched
S <- Xdfcs |> dplyr::select(prop_spring, prop_summer, prop_fall)

#A vector of the max distance searched from turbines in each report
dist_searched_transform <- as.numeric(Xdfcs$dist_searched_transform)

4. Prepare the model matrix and other model inputs

Using the centered and scaled covariate values, we set up the model.matrix of covariates for the model by selecting those to include. We also select the species to include in the model and subset the species ranges based on the selected species.

#Use an offset for the number of turbines relevant to each estimate
offset <- Xdf$n_turbine

#Set up the model matrix by selecting covariates to include in the analysis
X <- model.matrix(
  ~ 1 +
  mean_rsa +
  age +
  physdiv_av +
  mean_lrs +
  lrs_mean_windspeed_av +
  nlcd_ilr_1 +
  nlcd_ilr_2 +
  nlcd_ilr_3 +
  nlcd_ilr_4 +
  nlcd_ilr_5 +
  nlcd_ilr_6 +
  nlcd_ilr_7 +
  nlcd_ilr_8 +
  nlcd_ilr_9,
  Xdfcs)

#Remove the intercept from the model matrix, because a new intercept coefficient
#will be estimated for each species.
if(ncol(X) == 2){
  X <- as.matrix(X[,2], nrow = nrow(public_pcm_with_spc), ncol = 1)
} else {
  X <- X[,2:ncol(X)]
}

#Select the species to include in the model
species_columns <- c(
  "LACI (Hoary)",
  "LABO (Eastern red)",
  "EPFU (Big brown)"
)

#Subset the species range columns based on the selected species
R <- public_pcm_with_spc[,substr(species_columns, 1, 4)] |>
  st_drop_geometry()

Run the model

The data objects that were created above can be used as inputs to the multi-species model that we described in Model Details. The function-call below is set to a low number of iterations in order to demonstrate how the model runs in a quick fashion. However, in order to achieve convergence, a much higher number of iterations is suggested. In the comments, we provide a suggested number of iterations that we have used with this model to reach approximate convergence. Running the code below should take approximately 10 minutes on a standard 2020s computer. Expect the model to require up to several hours when increasing the number of iterations to achieve convergence.

Please keep in mind when running this model and reviewing the analysis, that the model is run only on the publicly available data included in this release. In the near future we plan to publish our results on the full data set, which includes business confidential reports. Our analyses have shown that the results on the full data set are more robust and don’t necessarily agree fully with the example results described below.

set.seed(20250404)

#Begin the model run with the processed inputs
#We recommend running the model for more iterations than specified in the demo
#Consider replacing the jags_pars with the following parameters
#These parameters will produce better convergence but require
#more runtime (several hours)
#   jags_pars = list(
#     ni = 1000000,
#     nc = 3,
#     na = 50000,
#     nb = 50000,
#     nt = 10,
#     parallel = TRUE),
system.time({
  jl <- NationalBatsAndWind::run_multispecies_model(
    public_pcm_with_spc,
    species_columns,
    X, 
    R,
    S,
    offset,
    adj_distance,
    dist_searched_transform,
    cs_save,
    jags_pars = list(
      ni = 100000,
      nc = 3,
      na = 50000,
      nb = 50000,
      nt = 3,
      parallel=TRUE),
    ecor = FALSE,
    model = "nmixture",
    family = "poisson",
    params = c(
      "beta_range", "beta_season", "beta_dist", "beta", "beta_mu",
      "eps", "eps_sigma"
    ))})

#Save the JAGS outputs in two objects
jags_data <- jl$jags_data
jags_fit <- jl$jags_fit

These outputs can also be saved to a file to load at another time.

#If desired, save the outputs to a file
save(jags_data, file = "../data/models/model_run_data.rda")
save(jags_fit, file = "../data/models/model_run_fit.rda")

Evaluate and visualize model outputs/predictions

Previous model outputs can be reloaded. If you are returning to this vignette, please run all code chunks prior to fitting and saving the model (the previous two code chunks), so that all R objects prior to the model run are available.

#If desired, load previous model outputs without having to rerun the models
#Make sure to run the other code chunks before the "Run the model" section also
load("../data/models/model_run_data.rda")
load("../data/models/model_run_fit.rda")

Evaluate convergence

Traceplots and whiskerplots can help evaluate convergence on model parameters. For a model run with a low number of iterations and a high number of parameters, it isn’t likely that convergence will be reached. Traceplots for each parameter should look like a “field of grass” that is centered around one mean across all iterations. Another sign of convergence is R-hat values < 1.1 or ideally < 1.01.

#Consider checking the traceplots of the following parameters for convergence.
#Edit the parameters variable in the traceplot() function below:
#parameters="beta_range"
#parameters="beta_dist"
#parameters="beta_season"
#parameters="beta"
jagsUI::traceplot(jags_fit, parameters = "beta_range")

Plot species intercept coefficient estimates

Our multi-species model of bat fatalities at wind energy facilities estimates a variety of parameters that can be summed together to predict species-specific fatality rates. First, we estimate the average annual rate of bats killed per turbine for each species. In the model, this parameter is \(\beta_s\). This parameter can be interpreted as the number of bats killed per turbine after marginalizing across the mean of all other covariate effects.

#Use the simulated parameter values to determine the mean and lower/upper bounds
#of the 95% credible interval for each intercept parameter
range_intercepts_data <- data.frame(
  covariates = species_columns, 
  color = jags_fit$overlap0$beta_range,
  lo = apply(exp(jags_fit$sims.list$beta_range), 2, quantile, 0.025),
  mean = apply(exp(jags_fit$sims.list$beta_range), 2, mean),
  up = apply(exp(jags_fit$sims.list$beta_range), 2, quantile, 0.975))

#Create a pretty plot of the intercept coefficients
range_intercepts_plot <- ggplot(range_intercepts_data) +
  geom_errorbar(aes(x = covariates, ymin = lo, ymax = up)) +
  geom_point(aes(x = covariates, y =  mean)) +
  geom_hline(yintercept = 0) +
  ylim(min(c(range_intercepts_data$lo, 0)), max(range_intercepts_data$up)) +
  coord_flip() +
  theme_minimal() +
  #ggtitle("Average annual fatalities per turbine") +
  theme(plot.background = element_rect(fill = "white", color = "white"),
        legend.position = "bottomleft",
        text = element_text(size=16, colour="black"),
        axis.text=element_text(size=12, colour="black"),
        axis.title=element_text(size=12,face="bold")) +
    ylab("Average annual fatalities per turbine") +
    xlab("Species")

#Save the plot for later
ggsave(
  "plots/multispecies_all_range_intercepts.png",
  plot = range_intercepts_plot,
  width = 6,
  height = 6,
  units = c("in"),
  dpi = 300,
  create.dir = T
)

#Print out the plot below
range_intercepts_plot

Plot covariate effect hypermeans across all species

We can also estimate the average effect of each covariate in the model across all species. These parameters are considered hyper-means and correspond to the \(\mu_{\text{covariate}_j}\) parameters in the model.

#Use the simulated parameter values to determine the mean and lower/upper bounds
#of the 95% credible interval for each covariate effect hypermean
covariate_data <- data.frame(
    covariates = colnames(jags_data$X), 
    color = jags_fit$overlap0$beta_mu,
    lo = apply(jags_fit$sims.list$beta_mu[,], 2, quantile, 0.025),
    mean = apply(jags_fit$sims.list$beta_mu[,], 2, mean),
    up = apply(jags_fit$sims.list$beta_mu[,], 2, quantile, 0.975)) |>
    mutate(Covariate = case_when(
        covariates == "mean_rsa" ~ "Mean rotor swept area",
        covariates == "age" ~ "Years since built",
        covariates == "mean_lrs" ~ "Mean lower rotor sweep",
        covariates == "physdiv_av" ~ "Physiographic diversity",
        covariates == "lrs_mean_windspeed_av" ~ "Mean windspeed",
        covariates == "nlcd_ilr_1" ~ "NLCD Ratio 1",
        covariates == "nlcd_ilr_2" ~ "NLCD Ratio 2",
        covariates == "nlcd_ilr_3" ~ "NLCD Ratio 3",
        covariates == "nlcd_ilr_4" ~ "NLCD Ratio 4",
        covariates == "nlcd_ilr_5" ~ "NLCD Ratio 5",
        covariates == "nlcd_ilr_6" ~ "NLCD Ratio 6",
        covariates == "nlcd_ilr_7" ~ "NLCD Ratio 7",
        covariates == "nlcd_ilr_8" ~ "NLCD Ratio 8",
        covariates == "nlcd_ilr_9" ~ "NLCD Ratio 9"
    ))

#Set the order of the covariates to appear in the plot
order <- c("Mean rotor swept area",
           "Years since built",
           "Mean lower rotor sweep", 
           "Physiographic diversity",
           "Mean windspeed",
           "NLCD Ratio 1",
           "NLCD Ratio 2",
           "NLCD Ratio 3",
           "NLCD Ratio 4",
           "NLCD Ratio 5",
           "NLCD Ratio 6",
           "NLCD Ratio 7",
           "NLCD Ratio 8",
           "NLCD Ratio 9")

#Create a pretty plot of the covariate effect coefficients
covariate_plot <- ggplot(covariate_data |> filter(Covariate %in% order)) +
    geom_errorbar(aes(x = Covariate, ymin = lo, ymax = up, color = color)) +
    scale_color_manual(values=c("FALSE"="red","TRUE"="blue")) +
    geom_point(aes(x = Covariate, y =  mean)) +
    geom_hline(yintercept = 0) +
    scale_x_discrete(limits = rev(order)) + 
    coord_flip() +
    ylim(min(covariate_data$lo), max(covariate_data$up)) +
    labs(color = "95% credible interval\noverlaps 0") +
    theme_minimal() +
    ggtitle("Effects over all species") +
    theme(title = element_text(size = 12, colour = "black", ),
          plot.background = element_rect(fill = "white", color = "white"),
          legend.position = "bottomleft",
          text = element_text(size=16, colour="black"),
          axis.text=element_text(size=12, colour="black"),
          axis.title=element_text(size=12,face="bold")) +
    ylab("Standardized effect size") +
    xlab("Covariate")

#Save the plot for later
ggsave(
  paste0("plots/covariate_effects_all_species.png"),
  plot = covariate_plot,
  width = 6,
  height = 6,
  units = c("in"),
  dpi = 300,
  create.dir = T
)

#Print out the plot below
covariate_plot

The hyper-means do not reveal much information. However, in the next section, we will consider how these parameters could have strong species-level effects.

Render species-specific predictions and visualizations

Create species-specific covariate effect plots

Let’s recreate the plot of covariate effects above, but plot how the effects vary for each species. These parameters are the species-specific effects in the model \(\beta_{\text{covariate}_j}\).

#Determine the species name (four letter code) from the species_columns
species_names <- substr(species_columns, 1, 4)

#For each selected species, create the plot
for(species in seq(1,length(species_names))) {

  #Select the species
  species_name <- species_names[species]
  
  #Determine the range of the species
  species_range <- spc_ranges |> filter(Spp_code == species_name)
  
  #Use the simulated parameter values to determine the mean and bounds
  #of the 95% credible interval for each species-specific covariate effect
  covariate_data <- data.frame(
    covariates = colnames(jags_data$X), 
    color = jags_fit$overlap0$beta[,species],
    lo = apply(jags_fit$sims.list$beta[,,species], 2, quantile, 0.025),
    mean = apply(jags_fit$sims.list$beta[,,species], 2, mean),
    up = apply(jags_fit$sims.list$beta[,,species], 2, quantile, 0.975)) |>
    mutate(Covariate = case_when(
        covariates == "mean_rsa" ~ "Mean rotor swept area",
        covariates == "age" ~ "Years since built",
        covariates == "mean_lrs" ~ "Mean lower rotor sweep",
        covariates == "physdiv_av" ~ "Physiographic diversity",
        covariates == "lrs_mean_windspeed_av" ~ "Mean windspeed",
        covariates == "nlcd_ilr_1" ~ "NLCD Ratio 1",
        covariates == "nlcd_ilr_2" ~ "NLCD Ratio 2",
        covariates == "nlcd_ilr_3" ~ "NLCD Ratio 3",
        covariates == "nlcd_ilr_4" ~ "NLCD Ratio 4",
        covariates == "nlcd_ilr_5" ~ "NLCD Ratio 5",
        covariates == "nlcd_ilr_6" ~ "NLCD Ratio 6",
        covariates == "nlcd_ilr_7" ~ "NLCD Ratio 7",
        covariates == "nlcd_ilr_8" ~ "NLCD Ratio 8",
        covariates == "nlcd_ilr_9" ~ "NLCD Ratio 9"
    ))

#Set the order of the covariates to appear in the plot
order <- c("Mean rotor swept area",
           "Years since built",
           "Mean lower rotor sweep",
           "Physiographic diversity",
           "Mean windspeed",
           "NLCD Ratio 1",
           "NLCD Ratio 2",
           "NLCD Ratio 3",
           "NLCD Ratio 4",
           "NLCD Ratio 5",
           "NLCD Ratio 6",
           "NLCD Ratio 7",
           "NLCD Ratio 8",
           "NLCD Ratio 9")
  
  #Create a pretty plot of the covariate effect coefficients
  spc_covariate_plot <- ggplot(covariate_data |> filter(Covariate %in% order)) +
    geom_errorbar(aes(x = Covariate, ymin = lo, ymax = up, color = color)) +
    scale_color_manual(values=c("FALSE"="red","TRUE"="blue")) +
    geom_point(aes(x = Covariate, y =  mean)) +
    geom_hline(yintercept = 0) +
    scale_x_discrete(limits = rev(order)) + 
    coord_flip() +
    ylim(min(covariate_data$lo), max(covariate_data$up)) +
    labs(color = "95% credible interval\noverlaps 0") +
    theme_minimal() +
    ggtitle(paste0(species_name, "-specific effects")) +
    theme(title = element_text(size = 12, colour = "black", ),
          plot.background = element_rect(fill = "white", color = "white"),
          legend.position = "bottomleft",
          text = element_text(size=16, colour="black"),
          axis.text=element_text(size=12, colour="black"),
          axis.title=element_text(size=12,face="bold")) +
    ylab("Standardized effect size") +
    xlab("Covariate")
  
  #Save the plot for later
  ggsave(
    paste0("plots/", species_name, "_all_covariate_effects.png"),
    plot = spc_covariate_plot,
    width = 6,
    height = 6,
    units = c("in"),
    dpi = 300,
    create.dir = T
  )
  
}

The plots that we generate in the code above are saved to the vignettes/plots folder. For each species, the code generates a “????_all_covariate_effects.png” file, where the question marks are replaced with each four-letter species code. Let’s look at the species-specific covariate effects for hoary bat (LACI = Lasiurus cinereus):

The species-specific results for hoary bat reveal stronger covariate effects. As the age of the facility and turbine ground clearance (mean lower rotor sweep) increase, there is a lower rate of hoary bat fatalities observed. Keep in mind, these are correlative observations, not necessarily causal and may be confounded with other covariates not included in the model. Additionally, there are a few landcover effects that are strongly associated with increased or decreased hoary bat fatality rates. These landcover effects are difficult to interpret directly, however, because they are based on a isometric logratio transformation which incorporates multiple landcovers into each covariate. In addition, please recall that this model evaluates only the publicly available data included in this release. In the near future we plan to publish our results on the full data set, which includes business confidential reports that we cannot share as raw data. Our analyses have shown that the results on the full data set are more robust and don’t necessarily agree fully with the example results described below. Finally, there may be some stochasticity in model runs, especially with a low number of iterations, so the strength of the effects described here may differ slightly.

Check out the additional ????_all_covariate_effects.png

Estimate the total annual and cumulative fatalities of each species

Using the estimated covariate relationships, we can estimate total bat fatalities for each species across a retrospective prediction for all wind turbines that have been on the landscape from 2001 to 2024. Recall that we compiled all turbines and wind facilities from the U.S. Wind Turbine Database into the NationalBatsAndWind::turbines and NationalBatsAndWind::facilities data frames in this package. For each of these facilities, we calculate the covariates used in the model above, then use the covariate effects to predict the total number of fatalities for each species that occurred at specific facilities over each year.

#Determine the species name (four letter code) from the species_columns
species_names <- substr(species_columns, 1, 4)

#For each selected species, create the plot
for(species in seq(1,length(species_names))) {

  #Select the species
  species_name <- species_names[species]
  
  #Determine the range of the species
  species_range <- spc_ranges |> filter(Spp_code == species_name)
    
  #Retrieve the center/scale information for the covariates
  cs_save <- jags_data$cs_save
  
  #Create a data.frame of parameter values to run predictions on
  parameter_grid <- data.frame(intercept = rep(1, 24),
                               prop_spring = rep(1, 24),
                               prop_summer = rep(1, 24),
                               prop_fall = rep(1, 24),
                               cutin = rep(0, 24),
                               max_cutin_speed = rep(3, 24),
                               year = c(2001:2024))
  
  #Edit sf settings to avoid warnings/errors
  sf::sf_use_s2(FALSE)
  
  #Determine which facilities fall within the selected species range
  facilities_st_in_range <- st_join(facilities, species_range) |>
    filter(!is.na(Spp_code)) |>
    sf::st_drop_geometry() |>
    left_join(NationalBatsAndWind::nrel_grid, "fid")
  
  #Make of copy of the nrel_grid to save outputs to
  nrel_grid_out <- NationalBatsAndWind::nrel_grid
  
  #Print some messages to provide progress feedback to user
  print(paste0("Processing species: ", species_name))
  
  #For each row in the parameter grid
  for(jj in 1:nrow(parameter_grid)){
     
    #Filter to facilities in the species range to those that existed that year
    facilities_st_filtered <- facilities_st_in_range |>
      filter(start_year <= parameter_grid$year[jj] & end_year > parameter_grid$year[jj])

    #Create covariate data to predict fatalities at each of these facilities
    facilities_for_pred <- facilities_st_filtered |>
      mutate(
        year = parameter_grid$year[jj],
        age = parameter_grid$year[jj] - start_year,
        n_turbine = uswtdb_p_tnum,
        p_cap = interp_p_cap,
        sum_t_cap = interp_p_cap,
        total_rsa = uswtdb_p_tnum * 3.14 * (interp_rd / 2)^2,
        total_mw = interp_p_cap,
        years_since_wns = parameter_grid$year[jj] - wns_yoa_av,
        years_since_wns_nonneg = if_else(years_since_wns <= 0, 0, years_since_wns),
        intercept = parameter_grid[jj,1],
        prop_spring = parameter_grid[jj,2],
        prop_summer = parameter_grid[jj,3],
        prop_fall = parameter_grid[jj,4],
        mean_rsa = 3.14 * (interp_rd / 2)^2,
        mean_hh = interp_hh,
        mean_lrs = interp_lrs,
        vapor_pr_av_sq = NA,
        cutin = parameter_grid[jj,5],
        max_cutin_speed = parameter_grid[jj,6])
     
    #Center and scale the corresponding variables that are in the model
    for(i in which(rownames(jags_data$cs_save) %in% colnames(jags_data$X))){
      #Drop units on the column if present
      if(class(facilities_for_pred[[rownames(jags_data$cs_save)[i]]])=="units"){
        facilities_for_pred[[rownames(jags_data$cs_save)[i]]] <- units::drop_units(facilities_for_pred[[rownames(jags_data$cs_save)[i]]])
      }
      #Center and scale the columns according to matches in cs_save
      facilities_for_pred <- facilities_for_pred %>% 
        mutate(!!rownames(jags_data$cs_save)[i] := (!!sym(rownames(jags_data$cs_save)[i]) - jags_data$cs_save[i,1] )/jags_data$cs_save[i,2])
    }
    
    #After center and scaling calculate the quadratic covariates
    facilities_for_pred <- facilities_for_pred |>
      mutate(vapor_pr_av_sq = vapor_pr_av^2)

    #Predict fatalities for each facility in facilities_for_pred
    temp <- facilities_for_pred %>% 
      NationalBatsAndWind::multispecies_predict(
        jags_data,
        jags_fit,
        species = species,
        resimulate = F,
        fast = F)
     
    predno_q2.5 <- paste0('pred_q2.5_', parameter_grid$year[jj])
    predno_q50 <- paste0('pred_q50_', parameter_grid$year[jj])
    predno_mean <- paste0('pred_mean_', parameter_grid$year[jj])
    predno_q97.5 <- paste0('pred_q97.5_', parameter_grid$year[jj])
    
    #Calculate the 0.025, 0.5, and 0.975 quantiles and mean of the posteriors
    #across all fid cells included in the species range from the nrel_grid
    out <- data.frame(fid = temp$cells_modeled)
    out <- out %>%
      mutate(q2.5 = apply(temp$pred, 1, quantile, 0.025)) %>%
      mutate(q50 = apply(temp$pred, 1, quantile, 0.5)) %>%
      mutate(mean = apply(temp$pred, 1, mean)) %>%
      mutate(q97.5 = apply(temp$pred, 1, quantile, 0.975)) |>
      group_by(fid) |>
      #If there are multiple facilities in one fid, add these values together     
      summarize(!!predno_q2.5 := sum(q2.5),
                !!predno_q50 := sum(q50),
                !!predno_mean := sum(mean),
                !!predno_q97.5 := sum(q97.5))
    
    #Add the predictions to the copy of the nrel_grid
    nrel_grid_out <- nrel_grid_out |>
      left_join(out, by = 'fid')
  }
  
  #Calculate cumulative posterior medians for each state across all years
  state_totals <- data.frame(nrel_grid_out) |>
     group_by(state) |>
     summarize(across(starts_with("pred_q50"), \(x) sum(x, na.rm = TRUE))) |>
     mutate(state = tolower(state))
  
  #Calculate lower bound of 95% credible interval for each year
  year_totals_q2.5 <- data.frame(nrel_grid_out) |>
       summarize(across(starts_with("pred_q2.5"), \(x) sum(x, na.rm = TRUE)))
  #Calculate median for each year
  year_totals_q50 <- data.frame(nrel_grid_out) |>
       summarize(across(starts_with("pred_q50"), \(x) sum(x, na.rm = TRUE)))
  #Calculate upper bound of 95% credible interval for each year
  year_totals_q97.5 <- data.frame(nrel_grid_out) |>
       summarize(across(starts_with("pred_q97.5"), \(x) sum(x, na.rm = TRUE)))

  #Prepare a table of annual lower bounds, medians, and upper bounds
  totals <- cbind(
    c("lower", "median", "upper"),
    rbind(as.vector(year_totals_q2.5[1,]),
          as.vector(year_totals_q50[1,]),
          as.vector(year_totals_q97.5[1,])))
  colnames(totals) <- c("year", 2001:2024)
  totals <- t(totals)
  colnames(totals) <- totals[1,]
  totals <- totals[-1,]
  totals <- data.frame(totals)
  totals$year <- 2001:2024
  
  totals$lower <- as.numeric(totals$lower)
  totals$median <- as.numeric(totals$median)
  totals$upper <- as.numeric(totals$upper)
  totals$year <- as.numeric(totals$year)
  
  #Create an empty list of plots for storage
  plots <- list()
  
  #For each selected year create a plot of total fatalities per state
  for(year in c(2001,2005,2010,2015,2020,2024)) {
     
    filter_year <- year
    
    #Select the facilities that exist on the landscape during filter_year
    facilities_to_plot <- facilities |>
      filter(start_year <= filter_year & end_year > filter_year) |>
      st_transform(crs = 4326)

    #Create a pretty plot of the total number of bat fatalities per state
    #Include state outlines and points for each wind facility
    output <- states |> 
      left_join(state_totals, by = "state") |>
      ggplot() +
      geom_polygon(aes(
        x = x, 
        y = y,
        group = group,
        fill = !!sym(paste("pred_q50", year, sep = "_"))),
        color = "white",
        linewidth = 0.25)  +
      ggthemes::theme_map() +
      labs(fill = "Total bats killed") +
      ggtitle(paste("Year: ", year, sep = "")) +
      scale_fill_gradientn(
        colours=c("navyblue", "darkmagenta", "darkorange1", "gold"),
        limits = c(1,max(state_totals |> dplyr::select(starts_with("pred_q50")))*1.1),
        labels = scales::unit_format(unit = "", big.mark = ","),
        na.value="gray50") +
      ggthemes::theme_map() +
      labs(fill = "Bat fatalities") +
      theme_minimal() +
      theme(
        plot.background = element_rect(fill = "white", color = "white"),
        legend.position = "right",
        text = element_text(size=16, colour="black"),
        axis.text=element_text(size=12, colour="black"),
        axis.title=element_text(size=12,face="bold")) +
      geom_point(
        data = data.frame(cbind(x = c(-100), y = c(40), total = c(NA))),
        aes(x = x, y = y, size = total, shape = "Wind facility"),
        colour = "gray75") + 
      geom_sf(
        data = facilities_to_plot,
        aes(shape = "Location", color="white"),
        color = "white", shape = 1, size = 0.25, show.legend = "point") +
      guides(
        shape = guide_legend(
          "",
          override.aes = list(shape = 1, size = 2, color = "black"))) +
      ylab("") +
      xlab("") +
      coord_sf(crs = crs(NationalBatsAndWind::nrel_grid))

    #Add plot to list of plots to use with gifski
    plots[[length(plots)+1]] <- output

    #Save plot for later
    ggsave(
      paste0("plots/", species_name, "_state_fatalities_", year, ".png"),
      plot = output,
      width = 9,
      height = 4,
      units = c("in"),
      dpi = 300,
      create.dir = T
    )
     
  }
  
  #Combine plots into a .gif
  library(gifski)
  png_files <- list.files(
    "plots/",
    pattern = paste0(species_name, "_state_fatalities_.*png$"),
    full.names = TRUE)
  gifski(
    png_files,
    gif_file = paste0(
      "plots/",
      species_name,
      "_state_fatalities_animation.gif"),
    width = 2700,
    height = 1200,
    delay = 1)

  #Calculate the cumulative sum of each statistic across years
  totals_cumulative <- totals |>
    mutate(
      lower_cumul = cumsum(lower),
      median_cumul = cumsum(median),
      upper_cumul = cumsum(upper))
  
  #Output the statistics to a file
  write.csv(
    totals_cumulative,
    file = paste0("plots/", species_name, "_totals_cumulative.csv"))

  #Create a ribbon plot of total annual fatalities
  annual <- ggplot(totals, aes(x=year, y=median)) + 
     geom_point() + 
     geom_line() + 
     geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.1) +
     scale_y_continuous(
       limits = c(1,max(totals$upper)*1.1),
       labels = scales::unit_format(unit = "", big.mark = ","),
       na.value="gray50") +
     theme_minimal() +
     xlab("Year") +
     ylab("Annual bat fatalities") +
     theme(plot.background = element_rect(fill = "white", color = "white"),
           legend.position = "right",
           text = element_text(size=16, colour="black"),
           axis.text=element_text(size=12, colour="black"),
           axis.title=element_text(size=12,face="bold"))
  
  #Save the plot for later
  ggsave(
     paste0("plots/", species_name, "_annual_fatalities.png"),
     plot = annual,
     width = 9,
     height = 4,
     units = c("in"),
     dpi = 300,
     create.dir = T
   )
  
  #Create a plot of total cumulative fatalities over time
  cumulative <- ggplot(totals_cumulative, aes(x = year, y = median_cumul)) + 
    geom_point() + 
    geom_line() + 
    geom_ribbon(
      aes(ymin = lower_cumul, ymax = upper_cumul),
      linetype = 2, alpha = 0.1) +
    scale_y_continuous(
      limits = c(1, max(totals_cumulative$upper_cumul) * 1.1),
      labels = scales::unit_format(unit = "", big.mark = ","),
      na.value="gray50") +
    theme_minimal() +
    xlab("Year") +
    ylab("Cumulative bat fatalities") +
    theme(
      plot.background = element_rect(fill = "white", color = "white"),
      legend.position = "right",
      text = element_text(size=16, colour="black"),
      axis.text=element_text(size=12, colour="black"),
      axis.title=element_text(size=12,face="bold"))
  
  #Save plot for later
  ggsave(
    paste0("plots/", species_name, "_cumulative_fatalities.png"),
    plot = cumulative,
    width = 9,
    height = 4,
    units = c("in"),
    dpi = 300,
    create.dir = T
  )
   
}

Let’s again consider just the hoary bat (LACI). But consider reviewing the figures for the other species in the vignettes/plots/ folder.

First, we can plot the estimated number of hoary bat fatalities per state over time. The animation shows how annual fatalities may change over time. Our predictions become more refined with our full data set and more covariates.

The model predicts that total fatalities have reached several hundred thousand hoary bats annually in recent years (95% credible interval).

Cumulatively, this example version of the model estimates that several million hoary bats have been killed due to wind energy nationwide since 2001.

Running the code chucks above will produce additional plots for other bat species, including the big brown bat (EPFU = Eptesicus fuscus) and the Eastern red bat (LABO = Lasiurus borealis), which are stored in the vignettes/plots/ folder:

????_state_fatalities_animation.gif ????_annual_fatalities.png ????_cumulative_fatalities.png

Render relative risk map for each species

Finally, we can consider how relative risk for each species may vary across the U.S. by pairing our predictions with projections from the National Renewable Energy Laboratory (NREL) that describe how much wind energy could potentially be added to the landscape under different scenarios.

We predict the number of nationwide bat fatalities that could occur in association with varying scenarios of wind energy build-out using the code below:

#Determine the species name (four letter code) from the species_columns
species_names <- substr(species_columns, 1, 4)

#For each selected species, create relative risk maps of potential fatalities
for(species in seq(1,length(species_names))) {

  #Select the species
  species_name <- species_names[species]
  
  #Print some messages to provide progress feedback to user
  print(paste0("Processing species: ", species_name))
  
  #Determine the range of the species
  species_range <- spc_ranges |> filter(Spp_code == species_name)

  #Load the center/scale information for the covariates
  cs_save <- jags_data$cs_save
  
  #Create a data.frame of parameter values to run predictions on
  parameter_grid <- data.frame(intercept = 1,
                               prop_spring = 1,
                               prop_summer = 1,
                               prop_fall = 1,
                               cutin = 0,
                               max_cutin_speed = 3,
                               year = 2024)
  
  #Create a subset of facilities added after 2020 representing new turbines
  facilities_post_2020 <- facilities |>
    filter(start_year >= 2020)
  
  #Print a progress note
  print("Intersecting NREL grid with species range... this could take 5+ minutes!")
  
  #Determine which nrel_grid cells fall within species range
  nrel_grid_in_range <- NationalBatsAndWind::nrel_grid |> st_filter(species_range)

  #For each NREL scenario, create a plot of relative risk
  for(scenario in c("limited", "reference", "open")) {
    
    #Print another note
    print(paste0("Creating map for ", scenario, " scenario"))

    #Open the corresponding scenario grid
    scenario_grid <- read.csv(paste0(
      "../data-raw/grid/nrel_grid/",
      scenario,
      "_2030_moderate_115hh_170rd_supply-curve.csv")) |>
      dplyr::mutate(fid = sc_point_gid)
    
    #Determine covariate values for hypothetical facilities at each grid cell
    nrel_turbines_summary <- nrel_grid_in_range |>
      left_join(scenario_grid, by = "fid") |>
      mutate(
        year = parameter_grid$year[1],
        age = 0, #what if we put a new facility at this location?
        p_cap = capacity_mw,
        sum_t_cap = capacity_mw,
        mean_ttlh = hub_height + 170/2,
        mean_hh = hub_height,
        total_rsa = n_turbines * (170 / 2)^2 * 3.14,
        mean_rsa = (170 / 2)^2 * 3.14,
        n_turbine = n_turbines,
        total_mw = p_cap,
        years_since_wns = year - wns_yoa_av,
        years_since_wns_nonneg = if_else(years_since_wns < 0, 0, years_since_wns),
        mean_lrs = hub_height - 170/2,
        lrs_mean_windspeed_av = windspeed_10_m_av + (2/3) * (windspeed_40_m_av - windspeed_10_m_av),
        shear_ttlh_minus_lrs_av = windspeed_200_m_av - (windspeed_10_m_av + (2/3) * (windspeed_40_m_av - windspeed_10_m_av)),
        total_rsa_pct = total_rsa,
        shape_ind = median(facilities_post_2020$shape_ind),
        para = median(facilities_post_2020$para),
        mean_dist_to_edge = median(facilities_post_2020$mean_dist_to_edge),
        intercept = parameter_grid[1,1],
        prop_spring = parameter_grid[1,2],
        prop_summer = parameter_grid[1,3],
        prop_fall = parameter_grid[1,4],
        vapor_pr_av = vapor_pr_av,
        vapor_pr_av_sq = NA,
        cutin = parameter_grid[1,5],
        max_cutin_speed = parameter_grid[1,6],
        physdiv_av = physdiv_av)
    
    #Center and scale the covariates according to cs_save
    for(i in which(rownames(cs_save) %in% colnames(jags_data$X))){
      if(class(nrel_turbines_summary[[rownames(cs_save)[i]]])=="units"){
        nrel_turbines_summary[[rownames(cs_save)[i]]] <- units::drop_units(nrel_turbines_summary[[rownames(cs_save)[i]]])
      }
      nrel_turbines_summary <- nrel_turbines_summary %>% 
        mutate(!!rownames(cs_save)[i] :=  (!!sym(rownames(cs_save)[i]) - cs_save[i,1] )/cs_save[i,2])
    }
    
    #After center/scaling, calculate the quadratic effects for a few covariates
    nrel_turbines_summary <- nrel_turbines_summary |>
      mutate(vapor_pr_av_sq = vapor_pr_av^2)

    #Predict across the new cells
    temp <- nrel_turbines_summary %>%
      NationalBatsAndWind::multispecies_predict(
        jags_data,
        jags_fit,
        species = species,
        fast = T)
    
    #Calculate the 0.025, 0.5, and 0.975 quantiles and mean of the posteriors
    #across all fid cells included in the species range from the nrel_grid
    predno_q2.5 <- paste0('pred_q2.5_', parameter_grid$year[1])
    predno_q50 <- paste0('pred_q50_', parameter_grid$year[1])
    predno_q97.5 <- paste0('pred_q97.5_', parameter_grid$year[1])
    out <- data.frame(fid = temp$cells_modeled)
    out <- out %>%
      mutate(!!predno_q2.5 := apply(temp$pred, 1, quantile, 0.025, na.rm = TRUE)) %>%
      mutate(!!predno_q50 := apply(temp$pred, 1, quantile, 0.5, na.rm = TRUE)) %>%
      mutate(!!predno_q97.5 := apply(temp$pred, 1, quantile, 0.975, na.rm = TRUE))

    #Add predictions to a new copy of the nrel_grid
    nrel_grid_pred <- nrel_turbines_summary
    nrel_grid_pred <- nrel_grid_pred |>
      left_join(out, by = 'fid')
    
    #Create a map of the relative risk of potential fatalities for the species
    potential_fatalities <- ggplot() + 
      geom_sf(
        data = nrel_grid_pred,
        aes(fill = if_else(
          pred_q50_2024 > 5000,
          5000,
          if_else(pred_q50_2024 < 1, 1, pred_q50_2024))),
        color = "transparent") +
        scale_fill_gradientn(
          trans = "log",
          colors = gist_ncar(100),
          labels = scales::unit_format(unit = "", big.mark = ",", accuracy = NULL),
          na.value = "gray25",
          breaks = c(5,50,500,5000),
          limits = c(1,5000)) +
      geom_polygon(
        data = states,
        aes(x = x, y = y, group = group),
        color = "black",
        fill = "transparent",
        linewidth = 0.25) +
      guides(shape = guide_legend(
        "Fatality report",
        override.aes = list(shape=16, size = 2))) +
      guides(size = "none") + 
      ggthemes::theme_map() +
      labs(fill = paste0("Annual ", species_name, "\nfatality potential")) +
      labs(size = NULL) +
      theme_minimal() +
      theme(
        plot.background = element_rect(fill = "white", color = "white"),
        legend.position = "right",
        text = element_text(size=16, colour="black"),
        axis.text = element_text(size=12, colour="black"),
        axis.title = element_text(size=12,face="bold")) +
        ylab("") +
        xlab("") +
        coord_sf(crs = crs(NationalBatsAndWind::nrel_grid))
     
     #Save the map for later
     ggsave(
         paste0(
           "plots/",
           species_name, "_",
           scenario, "_potential_bat_fatalities.png"),
         plot = potential_fatalities,
         width = 9,
         height = 4,
         units = c("in"),
         dpi = 300,
         create.dir = T
     )
     
   }

}

Again, let’s consider only the results for hoary bats in this vignette.

Under the “limited” NREL scenario, wind energy build-out is restricted by a high number of social, ecological, and economic constraints. Here, we plot hoary bat fatality potential across the United States. Fatality potential is calculated by predicting the potential maximum number of hoary bat fatalities within a specific grid cell, here, under the “limited” NREL scenario. Grayed out cells indicate that wind energy is restricted under this scenario and cannot be built in this location.

Under the “reference” NREL scenario, wind energy build-out is restricted by a moderate number of social, ecological, and economic constraints. Here, we plot hoary bat fatality potential for each grid cell under the “reference” NREL scenario.

Under the “open” NREL scenario, wind energy build-out has the lowest number of social, ecological, and economic constraints. Finally, we plot hoary bat fatality potential under the maximum amount of wind energy that might be built within each grid cell according to the “open” NREL scenario.

These maps meant to provide a general idea of where hoary bats may be at higher risk of fatalities due to wind energy build-out across the U.S.

Please review the plots that we have created for the other species, big brown bat (EPFU) and Eastern red bat (LABO). Check the following figures in the vignettes/plots/ folder:

????_limited_potential_bat_fatalities.png ????_reference_potential_bat_fatalities.png ????_open_potential_bat_fatalities.png

Here is the map for the “open” NREL scenario and annual potential big brown bat (EPFU) fatalities:

Finally, here is the map for the “open” NREL scenario and annual potential Eastern red bat (LABO) fatalities:

Again, please recall that this model evaluates only the publicly available data included in this release. In the near future we plan to publish our results on the full data set, which includes business confidential reports that we cannot share as raw data. Our analyses have shown that the results on the full data set are more robust and don’t necessarily agree fully with the example results we describe here. Additionally, we explore more covariate relationships, which refine our spatial predictions.